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Abstract. Some known results for locating the roots of polynomials are extended to the case of 
matrix polynomials. In particular, a theorem by A.E. Pellet [Bulletin des Sciences Mathematiques, 
(2), vol 5 (1881), pp. 393-395], some results of D.A. Bini [Numer. Algorithms 13:179-200, 1996] based 
on the Newton polygon technique, and recent results of M. Akian, S. Gaubert and M. Sharify (see 
in particular [LNCIS, 389, Springer p. p. 291-303] and [M. Sharify, Ph.D. thesis, Ecole Polytechnique, 
ParisTech, 2011]). These extensions are applied for determining effective initial approximations 
for the numerical computation of the eigenvalues of matrix polynomials by means of simultaneous 
iterations, like the Ehrlich-Aberth method. Numerical experiments that show the computational 
advantage of these results are presented. 
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1. Introduction. Consider a square matrix polynomial A(x) = Yl^o-A-iX 1 , 
where Ai are m x m matrices with complex entries and assume that A(x) is reg- 
ular, i.e., a{x) := dct^4(a;) is not identically zero. We recall that the roots of a(x) 
coincide with the eigenvalues of the matrix polynomial A(x) that is, the complex val- 
ues A for which there exists a nonzero vector v such that A(X)v = 0. Computing the 
eigenvalues of a matrix polynomials, known as polynomial eigenvalue problem, has 
recently received much attention [15] . 

In this paper we extend to the case of matrix polynomials some known bounds 
valid for the moduli of the roots of scalar polynomials like the Pellet theorem [2TJ 
123] , the Newton polygon construction used in 3J, applied in [3] and implemented 



in the package MPSolve (http://en.wikipedia.org/wiki/MPSolve) for computing 



polynomial roots to any guaranteed precision. We also extend a recent result by S. 
Gaubert and M. Sharify [8l [22] who shed more light on why the Newton polygon 
technique is so effective. Our results improve some of the upper and lower bounds 
to the moduli of the eigenvalues of a matrix polynomial given by N. Higham and 
F. Tisseur in [T31 Lemma 3.1]. 

1.1. Motivation. In the design of numerical algorithms for the simultaneous 
approximation of the roots of a polynomial a(x) = X)"=o xl(li w ith complex coeffi- 
cients, it is crucial to have some effective criterion to select a good set of starting 
values. In fact, the performance of methods like the Ehrlich-Aberth iteration [TJ [7] 
or the Durand-Kerner [SJ Q3] algorithm, is strongly influenced by the choice of the 
initial approximations j3] . A standard approach, followed in |^ , is to consider n 
values uniformly placed along a circle of center and radius r, say r = 1. This choice 
is effective only if the roots have moduli of the same order of magnitude. If there are 
roots which differ much in modulus, then this policy is not convenient since the num- 
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ber of iterations needed to arrive at numerical convergence might become extremely 
large. 

In [3] a technique has been introduced, based on a theorem by A.E. Pellet [2"Tll2"l] . 
[TH] and on the computation of the Newton polygon, which allows one to strongly 
reduce the number of iterations of the Ehrlich-Aberth method. This technique has 
been applied in [1] and implemented in the package MPSolve ( |http : //en . wikipediaTj 
org/wiki/MPSolve). This package, which computes guaranteed approximations to 
any desired precision to all the roots of a given polynomial, takes great advantage from 
the Newton polygon construction and is one of the fastest software tools available for 
polynomial root-finding. 

Let us introduce the following notation to denote an annulus centered at the origin 
of the complex plane: 

A(a, b) :={i€C, a < \x\ < b], (1.1) 

where < a < b. 

The theorem by A.E. Pellet, integrated by the results of [24] and [3], states the 
following property. 

Theorem 1.1. Given the polynomial w(x) = Y^i=o' WiX * with Wo,w n ^ 0, the 
equation 

n 

w K x K = J2 ( L2 ) 

i— 0, i^K 

has one real positive solution to if k — 0, one real positive solution s n if k = n, and 
either or 2 positive real solutions s K <t K if < k < n. Moreover, any polynomial 
b( x ) = S™=o ^ iX% suc h that \bi\ — \wi\, for i = 0, . . . ,n, has no roots of modulus less 
than to, no roots of modulus greater than s n , and no roots in the inner part of the 
annulus A(s K ,t K ) if s K and t K exist. Furthermore, denoting = ho < hi < . . . < 
hp — n the values of k for which equation (1.2) has two real positive solutions, then 



the number of roots of b{x) in the closed annulus A(t} l ._ 1 , Sh t ) is exactly hi — hi-x, 
fori=l,... ,p. 

The bounds provided by the Pellet theorem are strict since there exist polynomials 
w(x) with roots of modulus Sh i _ 1 and t^. Moreover, there exists a converse version 
of this theorem, given by J.L. Walsh [2"3I ITg] . 

The Newton polygon technique, as used in [3], works this way. Given a(x) — 
S"=o aiX ' ' wnere a 0i&n 7^ 0, the upper part of the convex hull of the discrete set 
{(i, log \ai\), i = 0, . . . ,n} is computed. Denoting hi, i = 0, q the abscissas of the 
vertices such that fco = < k\ < ■ ■ ■ < k q — n, the radii 

r i = \a ki _ 1 /a ki \ 1 /^- k *-i\ i = l,...,q (1.3) 

are formed and m, — hi — approximations are chosen in the circle of center 
and radius ri for i = 1, . . . , q. The integer ni; is called multiplicity of the radius r^. 
Observe that the sequence ri is such that r\ < r2 < • • • < r q and that — log are the 
slopes of the segments forming the Newton polygon. Figure |1.1| shows the Newton 
polygon for the polynomial a(x) = x 5 + 10 6 x 4 + x 3 + j^x 2 + 10 3 x + 1. 

The effectiveness of this technique is explained in |3J where the following result is 
shown. 

Theorem 1.2. Under the assumptions of Theorem \l.l\ it holds p < q, 



{h , h%,..., hp} c {k ,ki,...,k q }. 
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Fig. 1.1. Newton's polygon for the polynomial a{x) = x 5 + 10 6 x 4 + x 3 + j^qX 2 + 10 3 x + 1. 

Moreover, {n, . . . , r q }n]s hj , t hj [= for j = 1, . . . ,p. 

Therefore, the approximations chosen in the previously described way lie in the 
union of the closed annuli A(t} li _ 11 s^), i — 1, •• • ,p which, according to the Pellet 
theorem contain all the roots of all the polynomials having coefficients with the same 
moduli of the coefficients of a(x). Moreover, the number of initial approximations 
chosen this way in each annulus coincides with the number of roots that the polynomial 
has in the same annulus. The advantage of this approach is that the computation 
of the Newton polygon and of the radii is almost inexpensive since it requires 
O(n) operations, while computing the roots Sh t and i/, 4 is costly since it requires the 
solution of several polynomial equations. Recently, A. Melman (16j has proposed a 
cheap algorithm for approximating s K and t K . 

In [3] it is observed that any vertex (fc^log la^J) of the Newton polygon satisfies 
the following property 

Uki<v k ,, v ki = it/s i+1 = n+i, i = o,...,q—i, , , 

u ki := maxj< fei \aj/a ki \ 1/{ - ki ~ : '\ v ki := minj> fci \aj/a ki \ 1/( - ki ~^ ■ 

In certain cases, the radii given by the Newton polygon provide approximations 
to the moduli of the roots which are better than the ones given by the Pellet theorem. 
The closeness of the radii to the moduli of the roots of a(x) , holds in particular when 
the radii differ relatively much from each other, or, equivalcntly, when the vertices of 
the Newton polygon are sufficiently sharp. 

This has been recently proved by S. Gaubert and M. Sharify [HI HI]- In fact, using 
the theory of tropical polynomials, it turns out that the values 7"j, for i — l,...,q 
coincide with the so called tropical roots of a(x), and that the values rrii — fc, — 
are the multiplicities of the tropical roots r,-, for i = 1 , . . . , q. This fact is used in 
[8| 122] to prove the following interesting result. 

Theorem 1.3. // r j - 1 /r j , /r j + 1 < 1/9, fori — l,...,p, then any polynomial 
b(x) having coefficients with the same modulus of the corresponding coefficients of a{x) 
has ki — hi— i roots in the annulus A(ri/3, 3ri). 

That is, if three consecutive radii rj_i,rj and r^+i are sufficiently relatively far 
from each other, then the roots of any polynomial having coefficients with the same 
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modulus of the corresponding coefficients of a(x) are relatively close to the circle of 
center and radius r$. This explains the good performances of the software MPSolve 
where only the sufficiently sharp vertices of the Newton polygon are considered for 
placing initial approximation to the roots. 

An attempt to extend the Newton polygon technique to matrix polynomials is 
performed in [§] by relying on tropical algebra. The idea consists in associating with 
a matrix polynomial A(x) = ^2™= AiX l the Newton polygon constructed from the 
scalar polynomial w(x) = Y^,i=o II^MI 2 ^- An application of the results of [5] yields 
a scaling technique wich is shown to improve the backward stability of computing 
the eigenvalues of A(x), particularly in situations where the data have various orders 
of magnitude. The same idea of relying on the Newton polygon constructed from 
w(x) is used in [5] in the context of solving the polynomial eigenvalue problem with 
a root-finding approach. 

Moreover, it is proved in [5] for the quadratic matrix polynomial and in [221 Chap- 
ter 4] for the general case, that under assumptions involving condition numbers, there 
is one group of "large" eigenvalues, which have a maximal order of magnitude, given 
by the largest tropical root. A similar result holds for a group of small eigenvalues. 
Recently it has been proved in [2] that the sequence of absolute values of the eigen- 
values of A(x) is majorized by a sequence of these tropical roots, r^s. This extends to 
the case of matrix polynomials some bounds obtained by Hadamard [12] . Ostrowski 
and Polya [THU^O) for the roots of scalar polynomials. An attempt to extend Pellet's 
theorem to matrix polynomial by relying on the Gerschgorin disks has been performed 
by A. Melman [17] . 

1.2. New results. In this paper we provide extension to matrix polynomials of 
Theorems |1.2| |1.3| and arrive at an effective tool for selecting initial approxima- 
tions to the eigenvalues of a matrix polynomial. This tool, coupled with the Ehrlich- 
Aberth iteration, provides a robust solver for the polynomial eigenvalue problem. A 
preliminary description of an implementation of this solver is given in [5] . 

The Pellet theorem is extended by considering the equations 



n 

i— 0, i^K 



valid for all the k such that det A K ^ 0, which have either 2 or no positive solutions 
for k = 1 , . . . , n — 1 and the same equations for n = and n = n which have one 
positive solution if det Aq ^ or det A n ^ 0, respectively. 

The Newton polygon technique is extended by considering separately either the set 
of polynomials J27=o WA^AiWx 1 f° r K sucn that det A K ^ 0, or the single polynomial 
S"=o II^MI 2 ^- The latter case is subjected to the condition that Aj, i = 0, . . . ,n are 
well conditioned matrices. 



Theorem 1.3 is extended to matrix polynomials such that A,- t = Qi, QiQ* — atl 
for i = 0, . . . , n where the constants 3 and 9 are replaced by slightly different values. 
For general polynomials, computational evidence shows that the bounds deteriorates 
when the condition number of coefficients increases. 

The paper is organized as follows. In Section [2] we report the extensions of The- 
orems |1.1| \\.2\ |1.3| Section [3] contains the proofs of these extension. Finally, Section 
[4] describes the results of the numerical experiments that confirm the effectiveness of 
our extensions. 
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2. The main extensions. Throughout the paper A* denotes the conjugate 
transpose of the matrix A, p(A) is the spectral radius and \\A\\ — p(A* A) 1 / 2 is the 
2- norm. We denote by i the imaginary unit such that i 2 = —1. 

Define the class V m ,n of all the m x m matrix polynomials A(x) = Y^i=o AiX 1 
with At £ C« ixm 5 satisfying the following properties: 

• A(x) is regular and has degree n, that is, a(x) = det A(x) is not identically 

zero and A n ^ 0; 
. A ± 0. 

The latter condition is no loss of generality. In fact, in this case we may just consider 
the polynomial A(x)/x K , where k is the smallest integer such that A K ^ 0. 
Define also the class Q m „ C V m n such that 



S m ,„=£V, A i =a i Q i , QtQi=I, <7i>0}. (2.1) 



The class Q m ,n is given by all the matrix polynomials whose nonzero coefficients have 
unit spectral condition number. The expression aQ provides a first step in extending 
the complex number pe to a matrix, where a plays the role of p and Q of e l9 . 

The following result provides a first extension of the Pellet theorem. 

Theorem 2.1. Let A(x) e V m , n . 

1. If < k < n is such that det A K ^ then the equation 

x*= ]T WA^A^W (2.2) 

2—0, 

has either no real positive solution or two real positive solutions s K <t K . 

2. In the latter case, the polynomial a(x) — det A{x) has no roots in the inner 
part of the annulus A(s K ,t K ), while it has mtt roots of modulus less than or 
equal to s K . 

3. If k = and det Ao ^= then ( |2.2[ ) has only one real positive root to, moreover, 
the polynomial a(x) has no root of modulus less than to. 



4- If K = n and det A n =^ 0, then ( 2.2 I has only one real positive solution s n and 



the polynomial a(x) has no roots of modulus greater than s n . 

A consequence of the above result is given by the following corollary which pro- 
vides a further extension of Theorem 11.11 

Corollary 2.2. Let h < hi < ... < h p b e the values of k such that det A K ^ 



and there exist positive real solution(s) of (2.2). Then 
1- thi_! < s hl , i = 1, ■ ■ ■ ,p; 

2. there are m{hi — /ii-i) roots of a{x) in the annulus A{th i _ 11 Sh i ); 

3. there are no roots of the polynomial a(x) in the inner part of the annulus 
A(si li ,tf li ), where i — 0,1,..., p and we assume that Sh = s Q := 0,tf lp = 
t n := oo. 

Observe that in the case where m = 1, i.e., the matrix polynomial is a scalar 
polynomial, Corollary |2.2| coincides with Theorem |1.1| Moreover, the bounds to the 
moduli of the eigenvalues of A(x) given in the above results are strict since there exist 
matrix polynomials, say, polynomials with coefficients Ai — ail, which attain these 
bounds. 



Theorem 2.1 improves [131 Lemma 3.1] where the upper and lower bounds to the 
moduli of the eigenvalues of a matrix polynomial are given by the positive solutions 
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of the polynomial equations 



n 

1 = £11^1111 VII**- 

The improvement comes from the simple observation that ||AB|| < ||A||||B|| for any 
pair of square matrices A, B. 

If det^o = 0, clearly the lower bound to on the modulus of the roots of a(x) 
stated in part 3 of Theorem |2.1| is missing. In fact, there is at least one eigenvalue 
of A[x) equal to zero. Similarly, if dct A n = there is no upper bound s n to the 



modulus of the roots of a(x) stated in part 4 of Theorem 2.1 In fact, in this case 
there exist infinite eigenvalues. 

Notice that in Corollary |2.2| the value ho = exists if det Aq 0, and the 
value hp = n exists if det^4„ ^ 0. However, if the remaining coefficients A K , for 
K = 1, . ..,71— 1, even though non-singular, are very ill-conditioned, then it may 
happen that the set {hi, . . . ,/i p _i} is empty so that Corollary |2.2| does not provide 
much information. 

If A(x) G Qm.n, then the following extension of Theorem 1 1 . 1 1 holds . 

Theorem 2.3. Let A(x) € Q m<n so that Ai = o-j Qj and Q*Q% — I. Let Sh t ,t ht , 
i = l,...,p be the quantities given by Theorem 1.1 applied to w(x) — X)"=o aiXl • 



Then any matrix polynomial B{x) — 53t=o ^i^iX 1 G Qm.n 

1. m(hi — hi—x) eigenvalues in the annulus A{th i _ 1 , s^); 

2. no eigenvalues with modulus in the inner part of the annulus A(sh i ,th i )- 

2.1. The Newton polygon technique. The results given in the previous sub- 
section provide a useful tool for selecting initial approximations to the eigenvalues 
of A(x) to be refined by a polynomial root-finder based on simultaneous iterations. 
However, we may avoid to compute roots of polynomials and rely on the Newton 
polygon construction. 

In this section we provide some new results by using the Newton polygon tech- 
nique. We start by stating the following theorem. 



Theorem 2.4. Let A{x) € V m ,n- If k is such that (2.2) has two real positive 
solutions s K < t K , then u K < s K < t K < v K where 

u K ^max^JIA- 1 ^!! 1 /^), 

v K :=mnw ||Ar 1 Ai|l 1 /(«-*) 



If det Aq ^ then for k = (2.2 1 has a solution to and 



to < v :=mm\\A 1 A t \\- 1/l 

i>0 



If det A n ^ then for k = n (2.2) has a solution s n and 



>u n -maxHA- 1 ^!! 1 ^' 



Observe that for scalar polynomials the values u^. and w/ Ci are such that — 
Uk i+1 = ri, where are defined in (1.3) and ki are the abscissas of the vertices of the 
Newton polygon. 



The following result extends Theorem |1.2| to matrix polynomials. 
Theorem 2.5. Given A{x) e V m ,n, let S = : i = 0, —,q} be such that 
u K < Vk if and only if n G S , where u K and v K are defined in Theorem \2.J\ Then, 

1. p <q and {h , h p } C {fc , k q }; 

2. {{u 1 ,...,u q }U{v 1 ,...,v q }) n [s ki ,t ki ] = 0; 

3. v ki < u kt+1 , i = 0,...,q-l; 

4- if A(x) G Qm.nj then v ki — and v ki coincide with the vertices of the 

Newton polygon of the polynomial w(x) = X)"=o II II % l — X)"=o °~ xl ■ 



Therefore, the strategy of choosing m(fcj — ki-i) approximations placed along 
the circle of center and radius either rj = v ki or rj = u ki+1 is effective. In fact, 
these approximations lie in the union of the closed annuli A% of radii t) li l and Sft 4 , 
i = 1, . . . ,p, in the complex plane which, according to the extension of the Pellet 
theorem, contain all the eigenvalues of A(x). The computation of the radii Ti is 
cheap since it is reduced to compute the values fcj, i = 0, ... ,q defined in Theorem 



2.5 by evaluating the quantities u k . and v k . defined in Theorem 2.4 In the case of 



polynomials in Q m . n this computation is even cheaper since it is reduced to computing 
the Newton polygon of the polynomial w(x). 

Observe that for general matrix polynomials, u ki+1 does not generally coincide 
with v ki nor with the values obtained by computing the Newton polygon of w(x). 
However, in the practice of the computations, when the matrix coefficients corre- 
sponding to the vertices are well conditioned, there is not much difference between 
the values obtained in these different ways. 

The effectiveness of this strategy of selecting starting approximations is strength- 
ened by the following result which generalizes Theorem |1.3[ 

Theorem 2.6. Let A(x) = ^2 i=0 o~iQiX l € Q m ,n be a matrix polynomial of 
degree n and let ri, . . . , r g denote the radii of the Newton polygon associated with the 
polynomial w(x) = X)"=o °~i xl ■ Also, let ko, . . . , k q be the abscissas of the vertices of 
the Newton polygon and set mi = ki — There exist constants f,g such that 

12.11 < f < 12.12, 4.371 < g < 4.372 and 

1. for 1 < i < q, i/rj_i/r,-, ri/r i+ i < 1//, then A(x) has exactly mmi eigenval- 
ues in the annulus A(rij g,rig); 

2. for i — 1, if ri/ri+\ < 1/f then, A(x) has exactly mmi eigenvalues in the 
annulus A(ri/g' , rig), where g 1 = 2 + y/2; 

3. for i = q, if ri-i/ri < 1/f then, A(x) has exactly mmi eigenvalues in the 
annulus Afa/g, rig 1 ) . 



Observe that in the scalar case the values of the constants / and g are given by 
/ = 9 and g = 3 which are slightly better. 

3. The proofs. The key tool on which the proofs of our results rely is the 
generalization of Rouche theorem to the case of matrix polynomials provided in |10j , 
see also [53]. In this statement and throughout, we use the notation H y if the 
Hermitian matrix H is positive definite. 

Theorem 3.1. Let P{x) and Q(x) be square matrix polynomials, and T be a 
simple closed Jordan curve. If P{x)*P(x) — Q(x)*Q(x) y for all x £ T, then the 
polynomials p(x) := det(P(x)) and f(x) :— det(P(cc) +Q(x)) have the same number 
of roots in the open set bounded by T. 

We provide the proofs of the results listed in Section [2j Assume that det A K ^ 
and consider the matrix polynomial A(x) = Y^i=o Aix % which has the same eigen- 
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values as A(x). We start simply by applying Theorem 3.1 to the matrix polynomials 
P(x) = x K I and Q(x) — A(x) — P(x) where L is the circle of center and radius r. 

The condition P{x)*P{x) - Q{x)*Q{x) y turns into \x\ 2k I - Q{x)*Q{x) y 0. 
Moreover, since 

p(Q(x)*Q(x)) =\\Q(x)*Q(x)\\ < \\Q(x)*\\ ■ \\Q(x)\\ = \\Q(x)\\ 2 

NT, 



<( E \\ A * lA * 

4 — , i^K 



we deduce that the condition P(x)*P(x) — Q{x)*Q(x) >- is implied by 



\x\"> E \\A^A t 

i=0, i^K 

Thus, we may conclude with the following 



(3.1) 



Lemma 3.2. If (3.1 1 is satisfied for \x\ = r, then A(x) has run eigenvalues in 
the disk of center and radius r. 



Proof. If (3.1) is satisfied for |a;| = r, then P(x)P*(x) y Q(x)Q*(x) for P(x) 



3.1 



the matrix poly- 



x K I, Q{x) = A{x) — P(x) and \x\ = r. Therefore, by Theorem 
nomial P(x) + Q(x) = A(x) has as many eigenvalues of modulus less than r as the 
matrix polynomial P(x) = x K I, that is mn. □ 



We recall this known result which in [3] is proved by induction on k. 

Lemma 3.3. Let w(x) = YH=o w i xl > w o^ w n > 0, W{ > 0. The equation w K x K — 
S"=o i^K w i xl has only one real positive solution if k £ {0, n}, and either 2 or no real 
positive solutions if < n < n. 

ElLo.i^ Pk 1 ^!! -x*, in view of 



Applying Lemma 



3.3 



to the equation x K 
heorem 12.11 



Lemma [3721 one obtains 1 

Now, consider the set of indices W = {h < hi < 
and the equation 



< hp} such that det A K ^ 



E i 



A„, Ai 



X . 



has real positive solutions for k G %. Denote < t^, i = 0, . . . ,p these solutions, 
where we have set Sh = and th — oo. Observe that if dot Aq ^ then ho — and 
if det A n then h p — n. 



By applying Theorem |2.1| one deduces that the closed disk of center and radius 
Shi contains exactly mhi eigenvalues of A(x) for i — 1, . .. ,p, while there are none 
in the inner part of the annulus A{sh il th i )- This implies that Sh i+1 > tf ti% that is, 
part 1 of Corollary 2.2, and that there are m(hi — eigenvalues in the annulus 

A[thi- X i Shi)> i- e -; P ar ^ 2 of Corollary 2.2 Part 3 follows from a direct application of 



Theorem 2.1 so that the proof of Corollary |2.2| is complete. 

In the case where the matrix polynomial A(x) belongs to Q m ,m 
A~ x Ai = (o"i/<j K )J so that condition (3.1| turns into 



we find that 



xr > 



E 

2—0, i^K 



This way, the proof of Theorem 2.3 follows from Theorem 
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2.1 



3.1. Proofs related to the Newton polygon. Assume that < n < n and 



that equation (2.2) has real positive solutions s K < t K . Then, for any x € {s K , t K } one 
has 

n 

x K = WA-'A^x^WA-^Wx* 

2—0, i^K 

for any j ^ n. This implies that x > Tciaxj <K \\A~ 1 Aj\\ 1 / l - R ~^ =_u K and x < 

The cases 



mm 



A 1 j4 ? || 1 /( k ft — v K . This proves the first part of Theorem 



3>K W-^K 



2.4 



k = and n = n are treated similarly. 

Now consider Theorem |2.5| Parts 1 and 2 follow from Theorem |2.4| Concerning 
the inequality v ki < Uk i+ n we rely on the property \\H\\ > valid for any non- 

singular matrix H. In fact, for the sake of simplicity, denote k := ki and h := ki + \ so 
that k < h. Then 

u h =nmx\\A- 1 A J \\ 1 /^ > p^H 1 ^-*) > IKV 1 *) -1 !!" 1 ^"*' 



j<h 
n 

j>k " ' J j>k 



> minlKAr 1 ^)- 1 !!- 1 ^'-^ =min||^ 1 A i || 1 ^ fe -^ = v k . 



Concerning part 4, if A{x) € Q m .n, then in view of Theorem 2.3 the values of ki 
are the abscissas of the vertices of the Newton polygon for the polynomial w(x) = 
S"=o ° iX therefore v ki = Uk i+ i, and the proof of Theorem 2.5 is complete. 



3.2. Proof of Theorem 



2.6 



Let A(x) = J2i=o a iQi e Qm,n and w(x) = 
Y^i=o a i xt ■ Consider the Newton polygon corresponding to w(x) where ki-i,ki are 
the abscissas of two consecutive vertices corresponding to the ith edge. Also, let 
r±, . . . ,r q , denote the radii corresponding to different edges of the Newton polygon so 
that r\ < T2 < ■ ■ ■ < r q . Along the proof we refer to the radii as the tropical roots. 
Also, for the sake of notational simplicity we set k :— /c;_i and h :— ki. 

Let ri be the tropical root corresponding to the ith edge of the Newton polygon 
and consider the substitution y = rix. We define the matrix polynomial A(y) = 
Y^j=oAjyi as follows: 

n 

A(y) = (a k r*)- 1 C£ A i(riy) j ) ■ (3-2) 



Notice that A is an eigenvalue of A{y) if and only if r^A is an eigenvalue of A(x). The 
scaled matrix polynomial A(y) has the following property. 

Lemma 3.4 (Corollary of [2H Lemma 3.3.2]). Let A(y) be the matrix polynomial 
defined in (3.2 1 and let dj := \\Aj\\ for j = 0, . . . , n. Also let e :— ri-i/fi, 5 :— rijr^x, 



be the parameters measuring the separation between Ti and the previous and the next 
tropical roots, r,_i and rj+x, respectively. We have: 



e k-j jf j < / C; 

' 3 < \1 if k<j<h, 
S j - h if j > h. 



Proof. We only prove the first inequality since the other ones can be established 
by using a similar argument. Note that dj = {uk^i) -1 ^^. Due to [52} Lemma 3.3.2], 
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(7 j < Ufcr i _f for all < j < k. Thus, 



□ Now consider the decomposition of A(y) as the sum of two matrix polynomials: 



h k—1 n 

j = k j=0 j = h+l 



Remark 3.5. Notice that, although for the sake of notational simplicity we do 
not explicitly express this dependence, the definitions of A(y), P{y) and Q{y) depend 
on which edge of the Newton polygon we are considering. The following lemma 
provides the upper and lower bounds to the moduli of the eigenvalues of P(y). 

Lemma 3.6 (Corollary of [131 Lemma 4.1]). All the nonzero eigenvalues of P(y) 
lie in the annulus .4(1/2,2). 

Proof. Consider the polynomial y~ k P(y). Due to [HI Lemma 4.1] we have: 



(1 + Ili^H)- 1 min Hij-Vi < |A| < (1 + \\A^\\) max WA^-"^ 



where A is any eigenvalue of y k P(y). The result is established by applying Lemma 
|3.4| to the above inequalities. □ 

The idea of the proof is to look for the conditions on e = r^-i/r^ and S — rijr^x 
such that P(y)*P(y) — Q{y)*Q{y) >- holds on the boundaries of two disks of center 
zero and radi us r\ < 1/2 and r 2 > 2. Then, by the generalized Rouche theorem 



(Theorem 3.1 1, P{y) and A(y) will have the same number of eigenvalues inside these 



disks. Using Lemma 3.6 this implies that A(y) has m(h — k) eigenvalues which lie in 
the annulus A(ri, r 2 ); therefore, A{x) has m(h—k) eigenvalues which lie in the annulus 
Afari, TiV^)- This argument is akin to the one which is used in 22, Chapter 3] to 
prove Theorem |1.3[ valid for scalar polynomials. The proof relies on the following 
lemmas. 



Lemma 3.7. Let r :— \y\ and £ := h — k. Also, define the diagonal matrices 
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Di = dil m , i = 1, . . . , 3, and the Hermitian matrices Hi, i — 1 ... A, as follows: 



k-1 



DM := C^WUWm , D 2 (y) := H^H^X 

D 3 (y):=( X! 11^11 Wr» ; 

k<3i<32<h \k<ji<j2<h 



H 2 (y):= J2 A*A n {y*) n V n +\ E ^A(V*)V a 

0<il<j 2 <fe-l \0<ii<j 2 <fc-l 

H 3 (y)-.= E ^A(v*)V a +( E 44(!/7V z 

k— 1 n / fc— 1 n 

#4(y):=E E E E ^*A(y*)V 2 

Ji=Oj 2 =/i+l \ji=0j2=ft+l 

TTien we have: 

1. [P(y)]*P{y)-[Q{y)]*Q{y) = D 1 (y) + Hi(y)-D 2 (y)-D 3 (y)-H 2 (y)-H 3 (y)- 

Hi(y); 

2. ifr > 1, dt > r 2h r2 ~ 2+r " 



3. ifr< 1, di > • 
I \\Hx\\ < 2r k+1 > 



r 2 -l ' 
2k l-2r 2 +r 2e + 2 
1-r 2 



(r 2 -l)(r— 1) 

l-r5 



5. d 2 + d 3 + E- = 2 M < (e^F 1 + ^' 
Proof. The first equation is easily verified by direct computation. 
The proofs of inequalities 2. and 3. have been presented in }22l Lemma 3.3.5]. 
4. Note that 

||#i||<2 E Piillll^yi^ 1 ^ 2 < 2 E rjl+n (using Lemma\3^ 

k<ji<j2<h k<j!<j 2 <h 

h-1 

= 2(1 - r)- 1 E r 2jl+1 (l - r h -^) . 



Taking the sum over j\ completes the proof. 
5. It follows from Lemma HOI that 



k — 1 n 

d 2 <e 2fe E( r / e ) 2j ' ^3 < S- 2h E (<^) 2J > ll^ 2 ||<2e 2fc E (^A)^ 2 , 

o<ji<i 2 <fe-i 



k— 1 n 

32 



||ff 3 || < 26- 2h E W 1+ ' a , ll^4|| < 2e k 5- h E E ^AO* (»"<0 

h+l<j 1 <j 2 <n j 1 =0j 2 =h+l 
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Then we have 



fc-i 



d 2 + da + E M ^ I 6 " E^' + E 



i=2 



3=0 



j=h+l 



The proof is achieved by computing the right hand side of the above inequality. 
□ 

Lemma 3.8. For any edge i = 1, . . . ,q— 1 of the Newton polygon, if 8 — ri/r i+ i < 
1//, where 12.11 < / < 12.12, then A(y) and P(y) have the same number of eigenval- 
ues in the disk centered at zero and with radius g, where 4.371 < g < 4.372. For the 
last edge of the Newton polygon where i — q, A(y) and P(y) have the same number of 
eigenvalues in the disk centered at zero with radius g > 2 + y/2. 

Proof. Along the proof we assume that r > 1. Using Lemma |3.7[ a sufficient 
condition for [P{y)]*P{y) - [Q(y)]*Q{y) y is that 



di > d 2 + d 3 + p(-Hi + H 2 + H 3 + Hi 



Since p{-H x + H 2 + H 3 + H 4 ) < £* =1 \\Hi\\, (3.3 1 is implied by 



(3.3) 



d 1 -\\H 1 \\>d 2 



Eii^ 



(3.4) 



Observe that, by Lemma |3~7 



dx- \\HA >r 



2h i 



-4r + 2 2r- i+1 -r- 2 \ 



(r-1) 2 



(r-iy 



Thus, (3.4) is deduced from the following inequality 

-21 



4r + 2 2r- e+1 - 



> 



(r-1) 2 (r-1) 2 
Assume now that Sr < 1; since e < 1 we get 

1 - (Sr) 



Sr 



1 - (Sr) 



i—h 



r — e 



+ Sr 



1 - Sr 



< 



l-8r 



Sr 



Then (3.5) follows from 

r 2 - 4r + 2 2r 



-2f 



(r-1) 2 
which is equivalent to 
r 2 - 4r + 2 

2r 



(r-1) 2 



> 



r — 1 1 — Sr 



Sr 



r — 1 1 — Sr 



-21 ' 



-l _ i + _ 2( 5^+ 2 + £ r 



(r-1) 2 (r-l) 2 (l-5r) 
Now, assume that 8 < since ^ > 1, we have 
/+ 1 - 1 + - 2fr £ + 2 + Sr > r e+1 - 1 - 25r i+2 + 25r = (r 
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>( 



Sr 



1-Sr' 



(3.5) 



(3.6) 



(3.7) 



£+1 



l)(l-25r) > 



Then (3.7) follows from the following inequality 



4r 



> 



5 2 r 2 



(3.. 



(r-1) 2 (l-5r) 2 ' 

which can be written in the polynomial form 

p(6) := a{r)6 2 + b(r)S + c(r) > 0, 

where a(r) := r 2 (l — 2r), 6(r) := — 2rc(r), c(r) := r 2 — 4r + 2. Note first that r > 1 
implies a(r) < 0; let us consider the discriminant of p(S), i.e. A(r) := [&(r)] 2 — 
Aa(r)c(r) = 4r 2 (r — l) 2 c(r). We see that A(r) cannot be negative, otherwise p(<5) < 
Vr > 1. Thus, it must hold c(r) > 0, which implies that r > 2 + %/2- 
We deduce that p{5) > for all the values of 5 which satisfy 



6 < g+(r) := 



-c(r) + (r 



r(2r- 1) 



(3.9) 



The graph of <5+(r) is demonstrated in Figure 3.1 The maximum value of 6+(r) is 




6 18 30 



FlG. 3.1. The graph of <5+(r) as a function ofr. 



5max = Z± ^ 1 - \/18 



21-y/3 
2 



5 ¥ I + V 2 + Z # 

that [P(y)]*P(y) - 



0.08255 - 
4.3712. So for <5 < ^ max 



12.1136 



which is obtained at r = ro 



inequality (3 



holds, which implies 

[Q(y))*Q(y) >~ on t he b oundary of a disk of radius r ~ 4.3712. 
Then, by the Rouche theorem (Theorem 3.1 1, P{y) and A{y) = P(y) + Q{y) have the 



same number of zeros inside a disk of radius ro- This completes the proof of the first 
part of the lemma. 

Note now that for the last edge of the Newton polygon, since the terms D3, H3, H4 



are zero, inequality (3.6) becomes 



4r 



-2f 



(r-1) 2 



(r-1) 2 



> 



1 2 



r-1 



or equivalently, r 2 — 4r + 2 + 2r~ e+1 — 2r~ 2e > which holds for any r > 2 + y/2. □ 
Lemma 3.9. For any edge i = 2, . . . , q of the Newton polygon, if e — r^-i/rj < 
l/f, where 12.11 </< 12.12, then A{y) and P(y) have the same number of eigenval- 
ues in the disk centered at zero and with radius 1/g, where 4.371 < g < 4.372. For the 
first edge of the Newton polygon where i = 1, A(y) and P(y) have the same number 
of eigenvalues in the disk centered at zero, with radius 1/g where g > 2 + 
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Proof. Along the proof we assume that r < 1. We follow a similar argument to 
the previous lemma: the formulae that we will obtain are akin to those that we had 
gotten for the case r > 1. We will therefore give fewer details. 



Starting from (3.4) we get: 

„2Z+2 



2r A - 4r 



1 2r l+1 - r- 



(l-rf 



(1 



> 



l-jer- 1 ? ] 5rl+1 l_^(Sr) 



Sr 



(3.10) 



which is analogous to (3.5). Assume that e < r; we have 

2 



Sr 



i+i 



1 - (Sr) 



n—h 



1 - Sr 



< 



Sr 



i+i 



1 - Sr 



Thus, (3.10) is deduced from the following inequality: 
2r 2 



4r + 1 2r l+1 



r 2l+2 



(1-r) 2 
which is equivalent to 

2r 2 - 4r + 1 



(1- 



> 



< 



J+i 



-, 2 



■2r l+1 - 



(1-r) 2 

Now assume that e < £; we have 



r ; + 2 - 2e + er' +1 
(l-r) 2 (r-e) 



> 



r - r'+ 2 - 2e + er' +1 + er > r - r'+ 2 + 2er l+1 



(r- 



2e=(l-/ +1 )(r-2e)>0 



(3.11) 



(3.12) 



Thus, (3.12) is deduced from 



2r 2 - 4r + 1 



> 



(r- £ ) 2 



Setting p = r 1 we get ineq uality (3.8) and, following the arguments already used 

_ 7+3V3 



in the proof of Lemma 
0.08255 



18 



2 

0.22877 - 1 



we find a maximal value equal to e n 
— — r, which is obtained at r = r„ 1 = - V 3 + 2 '/ l 



disk of radius r. 



4.3712 

-1 



12.1136' io """""'^ " u 1 '0 2 

. For e < e max , [P(y)]*P(y) - [Q(y)]*Q(y) ^ on the boundary of a 



o 



Concerning the first edge of the Newton polygon, since the terms D2, H2, H4 are 



zero, (3.11) is replaced with 
2r 2 



4r + 1 2r l+1 



r 2l+2 



(1-r) 2 



(1- 



> 



-1 2 



1-r 



or equivalently, 2r 2 — 4r + 1 + 2r l+1 — 2r 2l+2 > 0, which holds for any r < 1 



V2 



□ 



Proof. [Proof of Theorem 



2.6 



3.6 



3.8 



3.9 



A(y) has m(/i — fc) 



By Lemmas 

eigenvalues lying in the annulus A(l/g,g), where 4.371 < g < 4.372. This fact implies 
that A(x) has m(h — k) eigenvalues in the annulus Airij g, gri), where rj denotes the 
tropical root corresponding to the ith edge of Newton polygon. □ 

Remark 3.10. The value of g that we obtained in Lemmas |3.8| and |3.9| yields 
uniform bounds, independent of the exact values of S and e. Yet, it is possible to 
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get tighter bounds if either S or e are smaller than the threshold value 1//. More 
precisely, due to inequality (|3.9[), the condition to be satisfied is 



-c(r) + (r - 1) y/rtrj - 5(2r 2 - r) > 

Suppose that one has a given value of 6, say S = 5q < S max = 1//. One can find the 
smallest r > 2 + v2 such that the above inequality holds. Notice that the function 



S+(r), defined as in (3.9), is increasing on the interval (2 + v2, ?"o]- Therefore, given 
any Sq < 1/ f there is a unique optimal radius r = f satisfying <5+(r) = 5q. 

As an example, when 6 — 0.05, the smallest r which satisfies the above inequality 
is r ~ 3.5142, while for S = 0.01 the smallest r is r ~ 3.4168 which is very close to 
2+V2. Following symmetric arguments, one one can show that when e is much smaller 
than the threshold value e max = 1/ /, the bound on the inner radius of the annulus 



can be improved. In this way, when either S or e are smaller than 1//, Theorem 2.6 
can be modified accordingly, providing sharper bounds for the eigenvalues of A(x). 

4. Numerical experiments. We have created a Matlab function which imple- 
ments the technique of choosing initial approximations to the eigenvalues of a matrix 
polynomial based on the results of Theorem |2.4| Even though this function is de- 
signed to deal with polynomials in the class Q m ,n, it can be applied to any matrix 
polynomial in V m .n- The function works in this way: the coefficient of the polynomial 
w(x) — J27=o are computed together with the associated Newton polygon 

which provides the values k t and r, = v^, i = l,...,q. Then m(k i+ i — kj) ini- 
tial approximations are uniformly placed in the circle of center and radius r^, for 
i = l,...,q. 

We have also implemented the Ehrlich Aberth iteration applied to the polynomial 
a(x) = det A(x) defined by 

(H-i) H N{x { f ) ) 
x, — x) , 1 = 1,..., mn, v = 0, 1, . . . , 

l-Mx^hV" l - 

1 "ft iL F i, 3 ^ >) >) 



N(x) = a(x)/a'(x) 



(4.1) 



starting from the initial approximations xf*\ i — 1, . . . , mn, where the Newton cor- 
rection N(x) is computed with the formula N(x) = l/tra.ce(A(x)~ 1 A'(x)). In this 
implementation, the iteration (4.1 1 is applied only to the components x " for which 



the numerical convergence has not occurred yet. We say that is numerically con- 
verged if either the Newton correction is relatively small, i.e., if |7V(a;^)| < e|a;^| for 
a suitable e > close to the machine precision, or if the reciprocal of the condition 
number of A(x^) is smaller than a given S > close to the machine precision. The 
execution is halted if either v — 5000 or if all the components xf have arrived at 
numerical convergence. 

Observe that each component xf' may converge with a number of iterations 



depending on i so that each simultaneous iteration in (4.1 1 does not have the same 



computational cost. In fact, while in the initial steps all the components in (4.1 ) must 



be updated, in the subsequent steps, when most of the components have arrived at 
convergence, only a few components xf ^ 1 must be updated. For this reason, it is not 
fair to compare performances by relying only on the maximum value reached by the 
parameter v which counts the number of simultaneous iterations. 
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m 


Newton polygon 
simul_it aver_it 


Unit circle 
simul_it aver_it 


5 


8 5.4 


243 


191 


10 


9 


5.5 


444 


375 


20 


11 


5.6 


855 


738 


40 


13 


6.1 


1594 


1466 



Table 4.1 



Number of simultaneous iterations and average iterations needed by the Ehrlich-Aberth 

iteration by choosing initial approximations on the unit circle or by using the Newton 

polygon technique. Polynomial with orthogonal coefficients scaled with constants (uj) = 
[1, 3.e5, 3.el0, l.el5, 0, 0, 0, 0, 0, l.e40, 0, 0, 0, 1] . 



Therefore, in our experiments besides counting the maximum number of simul- 
taneous iterations simul_it, that is, the maximum value reached by v, we have 
taken into account the average number of iteration per component aver_it, given 
by aver_it = ^2™^ Vi, where Vi is the number of steps needed for convergence 

of the ith component xf \ The quantity aver_it is more meaningful and represents 
the number of simultaneous iterations that one would obtain if all the components 
require the same number of iteration to arrive at convergence. The value of simul_it 
might be meaningful in a parallel environment where each processor can execute the 
iteration on a given component a;,. 

We have computed the values of simul_it and aver_it obtained by applying the 
Ehrlich-Aberth iteration starting with initial approximations xf^ uniformly placed 
along the unit circle and with initial approximations placed according to our criterion. 

The first set of experiments concerns matrix polynomials in the class Q m ,m i-e., 
polynomials with coefficients Ai — UiQi with Q*Qi — I and Oi > 0. The matrices 
Qi have been chosen as the orthogonal factors of the QR factorization of randomly 
generated matrices. Concerning the scalars at we have set 

cr=[l, 3.e5, 3.el0, l.el5, 0, 0, 0, 0, 0, l.e40, 0, 0, 0, 1] 

so that A{x) is a polynomial of degree 13 with eigenvalues of unbalanced moduli. We 
have chosen different values for m, more precisely, m = 5, 10, 20,40. Table |4~T] reports 
the number of iterations. It is interesting to point out that the reduction factor in 
the number of average and simultaneous iterations is quite large and grows almost 
linearly with the size m. 

We have also applied this technique to polynomials with randomly generated coef- 
ficients, which are not generally orthogonal, scaled by the same factors <Ji. The results 
are reported in Table [472] We may observe that the behavior is almost the same: the 
proposed strategy for choosing initial approximations still leads to a substantial de- 
crease of the number of iterations. It must be said that in the test performed, the 
condition numbers of the block coefficients corresponding to the vertices of the Newton 
polygon is not very large, the largest value encountered was around 5 . 0e3. 

5. Conclusions. Some known results valid for estimating the moduli of the roots 
of a polynomial have been extended to the case of matrix polynomials. These results 
have been applied to design a polynomial eigenvalue solver based on the Ehrlich- 
Aberth iteration. We have shown the effectiveness of our approach by means of 
numerical experiments. We plan to exploit these results in order to arrive at the im- 
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m 


Newton polygon 
simul_it aver_it 


Unit circle 
simul_it aver_it 


5 


9 


6.8 


240 


190 


10 


13 


7.7 


457 


372 


20 


16 


9 


851 


732 


40 


16 


10.4 


1597 


1457 



Table 4.2 



Number of simultaneous iterations and average iterations needed by the Ehrlich-Aberth 
iteration by choosing initial approximations on the unit circle or by using the New- 
ton polygon technique. Polynomial with random coefficients scaled by constants (cr;) = 
[1, 3.e5, 3.el0, 1.elB, 0, 0, 0, 0, 0, l.e40, 0,0, 0, 1]. 



plementation of a multiprecision matrix polynomial root-finder analogous to MPSolve 

a- 
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